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ABSTRACT 

Here we present an analysis of high-precision pulsar timing data taken as part of the North American 
Nanohertz Observatory for Gravitational waves (NANOGrav) project. We have observed 17 pulsars 
for a span of roughly five years using the Green Bank and Arecibo radio telescopes. We analyze these 
data using standard pulsar timing models, with the addition of time- variable dispersion measure and 
frequency-variable pulse shape terms. Sub-microsecond timing residuals are obtained in nearly all 
cases, and the best root-mean-square timing residuals in this set are ~30-50 ns. We present methods 
for analyzing post-fit timing residuals for the presence of a gravitational wave signal with a specified 
spectral shape. These properly and optimally take into account the timing fluctuation power removed 
by the model fit, and can be applied to either data from a single pulsar, or to a set of pulsars to 
detect a correlated signal. We apply these methods to our dataset to set an upper limit on the 
strength of the nHz-frequency stochastic supermassive black hole gravitational wave background of 
hc{l yr~^) < 7 X 10~^^ (95%). This result is dominated by the timing of the two best pulsars in the 
set, PSRs J1713+0747 and J1909-3744. 
Subject headings: pulsars, gravitational waves, . . . 



1. INTRODUCTION 

The direct detection of gravitational radiation (or grav- 
itational waves; GW) is currently a major goal in exper- 
imental physics. As described by general relativity, GW 
are freely propagating wave solutions to Einstein's equa- 
tion, or "ripples" in the space-time metric. Detecting 
GW would confirm another key prediction of general rel- 
ativity. GW are expected to be generated by nearly any 
configuration of accelerating mass, although due to the 
weakness of gravity, large masses or high accelerations 
are required to radiate significant GW. This means that 
astronomical objects are the only sources expected to 
produce measureable GW, and that we can in turn use 
these detections to learn about the GW sources them- 
selves; GW astronomy will provide an entirely new win- 
dow through which we can view the universe. Binary 
systems are expected to account for a large fraction of 
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the detectable GW signals, but we cannot discount "ex- 
otic" or unexpected sources either. 

One of the best tools we have for making precise as- 
tronomical measurements is the timing of radio pulses 
from spinning, magnetized neutron stars (pulsars). We 
receive a burst of radio waves once per spin period as the 
lighthouse-like beamed emission sweeps by Earth. The 
pulse times of arrival can be analyzed via a model that 
counts every single rotation of the star over years or even 
decades. This provides detailed information about the 
neutron star itself (via spindown rate and spin irregular- 
ities) , its binary orbit (via orbital Doppler and relativistic 
effects), and astrometry. Timing of double-neutron star 
binary systems has already provided strong evidence for 
the existence of GW, through measurements of orbital 
evolutio n induced by the generati on of GW by the binary 
system ([Tavlor fc Weisbeml 19891 ). The medium through 
which the radio pulses travel from the pulsar to Earth 
also affects the signal. This has been used extensively to 
probe the ionized compon ent of the inte rstellar medium 
in a variety of ways (e.g., IRickettiri990( ). The presence 
of GW along the line of sight also will affect the pulse 
travel times. This forms the basis for the use of pulsars 
as gravitational wave detectors. 

The infiuence of GW on pulsar timing, and its poten- 
tial use in GW detection, was fi rst noted over 30 years 
ago (jSazhinll 1 9 78t iDetweiledl 1 9 79f) . A majo r step forward 
was provided bv iHellings fc Downs! (|1983[ ). who showed 
that a GW signal will produce correlated timing fluctu- 
cations in a set of pulsars, a concept which came to be 
known as a pulsar timing array (PTA; iFoster fc Backeii 
119901) . A second major advance was the discovery of 
millisecond pulsars (|Backer et al.l 119821 ) - the high spin 
rate and stability of these objects improves GW sensi- 
tivity by orders of magnitude compared with canonical 
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'-^l-second pulsars. Pulsar timing is most sensitive to 
GW with period comparable to the observational time 
span, typically 1-10 years, so the GW frequency is in 
the nanohertz band. In this frequency range the bright- 
est expected sources are supermassive black hole binaries, 
whic h may be visible eithe r as a stochastic background 
(e.g., iJaffe fc BacAeil 120031: TS esana ct al. 2008) or as in- 
dividual sources (jSesana et al..,2009J . A numljer of pre- 
vious upper limits have been placed (iKaspi et al.lll995 
iJenet et all 12004 120061 : Ivan Haasteren et all 120111 ). but 
no successful detection has yet been accomplished. The 
topic continues to attract attention and a number of ded- 
icated pulsar timing array research projects are now ac- 
tive worldwide. 

In this paper, we present new high-precision timing 
observations of 17 millisecond pulsars, covering a five- 
year time span. These data were obtained as part of 
the North American Nanohertz Observatory for Grav- 
itational Waves (NANOGrav) proieciFI. using the two 
largest single-dish telescopes available, Arecibo Observa- 
tory and the NRAO Green Bank Telescope. In ^J2] we 
describe the observational setup and dataset. In ^J3] we 
present our methods and results for the determination 
of timing models. In 21 present a time-domain algo- 
rithm for measuring the cross-correlation between pairs 
of post-fit residuals, and apply this to our timing data 
to either detect or place limits on the stochastic gravi- 
tational wave background. We discuss the astrophysical 
implications of these results in fj5] and summarize our 
findings in Sj6l 

2. OBSERVATIONS 

The observations presented here were carried out over 
a five-year period, from 2005 to 2010, as part of a pro- 
gram specifically designed to measure or constrain the 
nHz- frequency stochastic gravitational wave background. 
The sources were observed with one or both of the 
305-m NAIC Arecibo Observatory or the 100-m NRAO 
Green Bank Telescope (GBT), with a typical observa- 
tional cadence of 4-6 weeks between sessions (see Fig.[T]). 
Scheduling at each observatory was independent, and the 
sessions were typically not performed simultaneously at 
the two telescopes. 

These sources were selected from the known population 
of bright millisecond pulsars (MSPs), excluding those 
found in globular clusters. Source selection was per- 
formed with the goal of obtaining the highest timing pre- 
cision possible, to maximize sensitivity to gravitational 
waves. Due to its higher gain, Arecibo was used to ob- 
serve all sources in its visible portion of the sky (^ 0° 
to 30° declination) , while sources outside this range were 
observed with the GBT, down to its minimum declina- 
tion of ~ —45°. One source presented in this paper, 
J1713-I-0747, was observed with both telescopes. [3 

Each pulsar was observed at two widely separated fre- 
quencies in order to track variation in dispersion measure 
(DM; see ^3.2\i . As the different frequency measurements 
use different receivers, these observations are also not si- 

http://www.nanograv.org 

1^ Another bright MSP, PSR B1937-I-21, is also observed with 
both telescopes as part of this project. We have for now excluded 
it from thi s analysis due to it s well-documented high levels of spin 
noise (e.g.,|K aspi et al."1994j) and interstellar medium systematics 
fe.g.. IRamchand ran et al. 200B ). 



multaneous. The choice of receivers used to observe a 
certain pulsar was determined from the availability and 
performance of the receiver systems at each telescope, 
and from the spectral index of the pulsar. At Arecibo, 
multi-frequency data are obtained in a single observing 
session, within ^1 hour for each source. At the GBT, the 
different frequencies were observed up to ~1 week apart. 
In a given epoch the observing time per frequency band 
for each source was anywhere from 15 to 45 minutes, with 
shorter integration times in general at Arecibo. The tele- 
scopes and frequencies used to observe each source, along 
with its basic physical parameters such as spin period and 
DM, are listed in Table[T] Three pulsars in this set (PSRs 
J1853-I-1308, J1910-hl256, and B1953-F29) were origi- 
nally observed as par t of a different project at Arecibo 
(jGonzalez et al.l[2011t ) and consequently only have single- 
band data for most of the time span. All receiver systems 
used at 800 MHz and above are sensitive to orthogonal 
linear polarizations. The 430 MHz system at Arecibo is 
sensitive to dual circular polarizations. 

All observations were performed using the identical As- 
tronomical Signal Processor (ASP) and Green Bank As- 
trono mical Signal Pro cessor (GASP) pulsar backend sys- 
tems (|Demorestl 120071 ) at Arecibo and the GBT respec- 
tively. These systems perform real-time coherent dedis- 
persion, full-Stokes detection, and pulse period folding in 
software, using a ~20-node Linux-based computer clus- 
ter. Channelized voltage data are supplied to the clus- 
ter from a SERENDIP-V FPGA boartO controlled by a 
Compact PCI single-board computer. The SERENDIP- 
V board performs analog-to-digital conversion and a dig- 
ital polyphase filterbank operation to split a 128-MHz 
bandwidth, dual-polarization signal input from the tele- 
scope into 32 4-MHz wide channels. Both the initial dig- 
itization of the signal and the channelized output use 
8-bit complex quantization. Further details of the hard- 
ware and real-time software used i n these observations 
were described bv iDemorestI ()2007D . While the 4-MHz 
channelization provides a minimum time resolution of 
250 ns, the final time resolution depends on the pulse pe- 
riod and number of profile bins used to average the profile 
- in this case, either 2048 (GBT) or 4096 (Arecibo) bins 
were used. Profiles were typically integrated for 1 minute 
(Arecibo) or 3 minutes (GBT). The total bandwidth pro- 
cessed varied with pulsar and observing frequency, lim- 
ited by either the real-time computational load or the 
receiver bandpass. A maximum of 64 MHz was used in 
most cases, with smaller bandwidths down to ^20 MHz 
used for low-frequency or high-DM observations. 

3. TIMING ANALYSIS 

In this section, we describe the procedures used in the 
timing portion of the data analysis. These can be split 
into two main areas: analysis dealing with pulse profiles, 
including polarization calibration and determination of 
pulse times of arrival fi 33.1|) : and fitting a physical timing 
model to the arrival times (i 33.2|) . The further analysis 
step of determining GW background limits is discussed 
later, in 21 

3.1. Calibration and Time of Arrival Estimation 
|https : //casper .berkeley ■ edu/galf a/"| 
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Fig. 1. — Overview of timing residuals for all sources, showing observational cadence and coverage during the five-year time span. The 
gap in 2007 was due to an extended maintenance period at both telescopes. The full scale of the y-axis is 10 /is in all cases. 
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TABLE 1 

List op observed millisecond pulsars: Basic parameters and observing setups. 



Source 


P 

(ms) 


dP/dt 
(10-20) 


DM 

(pc cm^'^) 


Pb 
(d) 


327 MHz 


Average flux density ( 
430 MHz 820 MHz 


mjy)° 
1.4 GHz 


2.3 GHz 


Obs 


J0030+0451 


4.87 


1.02 


4.33 






13.7 




1.4 




AO 


J0613-0200 


3.06 


0.96 


38.78 


1.2 






5.3 


2.0 




GBT 


J1012+5307 


5.26 


1.71 


9.02 


0.6 






7.6 


3.9 




GBT 


J1455-3330 


7.99 


2.43 


13.57 


76.2 






2.0 


1.1 




GBT 


J1600-3053 


3.60 


0.95 


52.33 


14.3 






3.1 


2.3 




GBT 


J1640+2224 


3.16 


0.28 


18.43 


175.5 




10.8 




1.0 




AO 


J1643-1224 


4.62 


1.85 


62.42 


147.0 






12.3 


4.2 




GBT 


J1713+0747 


4.57 


0.85 


15.99 


67.8 






8.8 


6.3 


3.6 


AO,GBT 


J1744-1134 


4.07 


0.89 


3.14 








7.6 


2.6 




GBT 


J1853+1308 


4.09 


0.87 


30.57 


115.7 








0.2 




AO 


B1855+09 


5.36 


1.78 


13.30 


12.3 




14.0 




4.0 




AO 


J1909-3744 


2.95 


1.40 


10.39 


1.5 






3.4 


1.4 




GBT 


J1910+1256 


4.98 


0.97 


34.48 


58.5 








0.2 




AO 


J1918-0642 


7.65 


2.57 


26.60 


10.9 






4.5 


1.8 




GBT 


B1953+29 


6.13 


2.97 


104.50 


117.3 








1.0 


0.1 


AO 


J2145-0750 


16.05 


2.98 


9.03 


6.8 






12.3 


3.2 




GBT 


J2317+1439 


3.45 


0.24 


21.90 


2.5 


32.2 


5.4 








AO 



" The presence of flux density values indicate the frequencies at which each source is observed. 
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As discussed in ^ the data products resulting from 
an observing session are a set of full- Stokes pulse pro- 
files integrated over 4 MHz radio bandwidth and 1-3 
minutes of time, into either 2048 or 4096 pulse phase 
bins. Following standard pulsar data analysis proce- 
dures, we aim to determine from the profile data a set 
of pulse times of arrival (TO As), i.e. times at which the 
apparent rotational phase of the pulsar passes through 
some fiducial point. This process involves several major 
steps: polarization calibration, template profile creation, 
additional profile averaging, and finally TOA measure- 
ment. For this work, we have performed all of these 
steps using two independent versions of pulsar data pro- 
cessmg s oftware, the first b ased on the PSRCHIVEEB 
software (iHotan et al.' '2004^ . and the second based on 
ASPFitsReader (|Ferdman 2008). This procedure pro- 
vides an important cross-check for errors in the analysis 
software that might otherwise be hard to detect. Both 
analysis pipelines performed the same procedures, aside 
from template determination which was done only via 
PSRCHIVE. 

Pulsar radio emission is typically highly polarized. 
While our current analysis relies only on the total in- 
tensity profiles for TOA determination, due to the high 
degree of polarization, calibration errors can still dis- 
tort t he intensity profile shape, leading to a TOA bias 
(e.g., Ivan StratenI 120061) . A complete description of 
the instrumental response to a polarized signal is pro- 
vided by the Mueller matrix, which is a radio- frequency- 
dependent linear transfor mation from the intrinsic to ob- 
serve d Stokes parameters (jHeiles et al.ll200ltlvan StratenI 
120041) . While we plan in future work to apply full Mueller 
matrix calibration to these data, for the current analysis 
we correct only the leading order terms, differential gain 
and phase between the two polarization components of 
the telescope signal. This is done via an injected cali- 
bration signal; immediately before or after each pulsar 
observation, a noise diode switched at 25 Hz is coupled 
into both polarization signal paths and measured with 
the pulsar backends. This provides a constant reference 
power versus frequency that is used to scale the pul- 
sar data. The equivalent flux density of the calibration 
signal is determined separately for each polarization by 
observing the noise diode along with a bright, unpolar- 
ized quasar of known flux density (B1442+101 at Green 
Bank, J1413-K1509 at Arecibo). This then provides a 
second scaling to convert the pulsar data to flux den- 
sity units (Jy) separately in each linear (or circular for 
certain receivers) polarization. The two calibrated total 
power terms are then added together to form the total 
intensity ("Stokes I") proflle, and the polarized profiles 
are not used further in this analysis. 

The calibrated pulse profiles are used to determine 
pulse TOAs and their uncertainties by fitting for a pulse 
phase shift bet ween each pr ofile and a standard "tem- 
plate" profile (|TavloHfl99a) . The template ideally is 
a noise-free representation of the average pulse proflle 
shape; any noise present in the template, especially if 
correlated wi th noise in the proflles, can bias the TOAs 
(jHotan et al.l 12005). For this work we determine tem- 
plate profiles from our measured data in a two-step 
process: the profiles are first roughly aligned using a 

|http : //psrchlve . sourcef orge . net | 



single-Gaussian template, then are summed together us- 
ing weights to optimize th e S/N in the final full-sum 
profile fsee iDemores^ 120071 Eqn. 2.10). We then ap- 
ply a tr anslation-invariant wavelet transform and thresh- 
olding ([Coifman &: Donohol I1995D to the profile to re- 
move its noise (Figure [5]) . The procedure is then iter- 
ated once, now using the new template for alignment, to 
produce the final template profiles used for determining 
TOAs. The wavelet noise removal procedure is imple- 
mented in the PSRCHIVE program psrsmooth. In this 
way, we obtained one template per pulsar per receiver 
used to observe it. While all other processing steps in 
this section were performed independently by both soft- 
ware pipelines, for consistency the same set of templates 
was used in both cases. 

The standard template-matching procedure used for 
TOA determination in this analysis is kno wn to suffer 
probl ems when applied to low-S/N data (jHotan et al.l 
I2005D . Therefore, it is useful to average profiles over 
as much time and radio bandwidth as possible to max- 
imize the S/N ratio before forming TOAs, while still 
retaining enough resolution to measure all appropriate 
instrumental or astrophysical effects. In this analysis 
we will retain the native instrumental frequency reso- 
lution (see tj3.2p . and have chosen to average over time 
all profiles in a given frequency channel from each ob- 
serving epoch (typically 30 minutes). TOAs are then 
measured for each average profile, resulting in a set of 
^-^20-30 TOAs (one per 4-MHz frequency channel) per 
each dual-receiver pair of observing epochs. The total 
number of TOAs for each pulsar are listed in Table [51 
As previously mentioned, we computed TOAs using two 
independent analysis pipelines. After verifying that the 
two pipelines produced consistent results, we focus the 
remainder of the analysis on the PSRCHIVE-produced 
data; all further results presented in this paper are spe- 
cific to these data. These TOAs are the inputs for the 
next part of the analysis procedure, fitting the timing 
model. The entire set of TOAs used in this analysis can 
be obtained as an electronic supplement to this paperF^ 

3.2. Timing Model Fit 

The second part of the timing analysis is to fit the 
measured pulse TOAs for each pulsar to a physical tim- 
ing model. The timing model predicts the apparent rota- 
tional phase of a pulsar based on a set of physical parame- 
ters describing the star's rotation (spin period, spin-down 
rate) , astrometry (position, proper motion, parallax) , bi- 
nary orbital motion, and general relativistic effects such 
as Shapiro delay. The model prediction is compared to 
the measured TOAs and best-fit parameter values are 
determined via minimization. This procedure is a 
fundamental part of pulsar astronomy and has been de- 
scribed many t imes in the literature (see for example 
iLorimer fc Kra mer 2004) . For all results presented here, 
we use the standard TEMPCEl timing analysis software. 
Preliminary comparis ons with analysis using the newer 
TEMP0^3 package (|Hobbs et al.ll2006l ) produce nearly 
identical results and will be presente d in future works 
(|Perrodin et al.ll20T2l:lEriis et al.ll20ll ). 

|http : //www. cv . nrao ■ edu/$\slm $pdemores/nanograv_data| 
'http : //tempo . sour cef orge . net 
http://tempo2.sourceforge.net 
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Fig. 2. — Full-sum profile and template profile for J1713-I-0747 at 1400 MHz. The top panel shows the full-sum profile at full scale. 
The middle panel shows the full-sum profile (points) and wavelet- denoised template version (line). The bottom panel shows the residual 
difference between the full-sum profile and template. 



TABLE 2 

Overview and results from timing model fits. 



Source 


#of 


# of parameters 




RMS 


Fit 


Epoch- 


averaged RMS {^^s)'^ 




TOAs" 


DM 


Profile Other'' 


(ms) 




Low-band'' 


High-band 


Combined 


J0030-I-0451 


545 


20 


26 


7 


0.604 


1.44 


0.019 


0.328 


0.148 


J0613-0200 


1113 


34 


45 


12 


0.781 


1.21 


0.021 


0.519 


0.178 


J1012-(-5307 


1678 


52 


53 


14 


1.327 


1.40 


0.192 


0.345 


0.276 


J1455-3330 


1100 


37 


53 


12 


4.010 


1.01 


0.363 


1.080 


0.787 


J1600-3053 


625 


21 


31 


14 


1.293 


1.45 


0.233 


0.141 


0.163 


J1640-I-2224 


631 


23 


26 


12 


0.562 


4.36 


0.057 


0.601 


0.409 


J1643-1224 


1266 


40 


48 


13 


2.892 


2.78 


0.589 


1.880 


1.467 


J1713-I-0747 


2368 


50 


111 


15 


0.106 


1.48 


0.092 


0.025 


0.030 


J1744-1134 


1617 


54 


49 


7 


0.617 


3.58 


0.139 


0.229 


0.198 


J1853-(-1308 


497 





34 


12 


1.028 


1.16 


0.271 


0.096 


0.255 


B1855-I-09 


702 


29 


21 


14 


0.395 


2.19 


0.277 


0.101 


0.111 


J1909-3744 


1001 


31 


37 


14 


0.181 


1.95 


0.011 


0.047 


0.038 


J1910-I-1256 


525 





34 


14 


1.394 


2.09 


0.712 


0.684 


0.708 


J1918-0642 


1306 


49 


37 


12 


1.271 


1.21 


0.129 


0.211 


0.203 


B1953-I-29 


208 





27 


12 


3.981 


0.98 


1.879 


0.543 


1.437 


J2 145-0750 


675 


20 


37 


12 


1.252 


1.97 


0.068 


0.494 


0.202 


J2317-I-1439 


458 


30 


12 


15 


0.496 


3.03 


0.373 


0.150 


0.251 



" One TOA per frequency channel per epoch. 

* "Other" parameters are all spin, astrometric and binary parameters as described in i|3.2l 
RMS computed from residuals averaged down to one point per receiver per epoch. See text for details. 
Note that in these results, the low-frequency RMS tends to be suppressed due to the DM(t) fit. 
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As the main goal of this analysis is to detect or limit 
the nHz gravitational wave background, a detailed dis- 
cussion of each pulsar's timing model parameters and the 
astrophysical significance of the results will not be pre- 
sented here. The TEMPO parameter files containing the 
model parameters and fit results can be obtained along 
with the TOAs in the electronic supplement to the pa- 
per. However, the overall strategy for fitting the timing 
models can be described as follows: 

1. The average pulsar spin frequency and frequency 
derivative (spin-down) were always fit parameters. 
No higher frequency derivatives were included. 

2. All five astrometric terms - two sky coordinates, 
two components of proper motion, and parallax - 
were always fitted parameters. 

3. Binary systems were fit for all five Keplerian 
parameters, using ether the "ELLl" or "DD" 
(jPamour fc DerueM I1985D timing models as ap- 
propriate. Additional relativistic or secular orbital 
terms were added only if the timing was signifi- 
cantly improved. 

In addition to the model parameters just described, 
the timing fit also included terms to correct for time- 
variable DM and frequency dependence of the pulse pro- 
file shape. Dispersion is a propagation effect caused by 
the travel of the radio pulses through the ionized inter- 
stellar medium (ISM). This causes a radio frequency- 
dependent shift in pulse arrival time proportional to 
DMiy~^. Astrophysically, DM represents the integrated 
column density of free electrons along the line of sight 
to the pulsar. Due to relative motion of the Earth and 
the pulsar, the effective path through the ISM changes 
with time, hence the DM vari es in a stochastic manner 
(e.g. iRamchandran et al.ll2006D . For this experiment, we 
obtained (non-simultaneous) dual-frequency data as de- 
scribed in f}2l specifically to measure and remove this ef- 
fect from the timing. We have done this by including a 
piecewise-constant DM(t) function in the same fit that 
determines the other timing model parameters. For each 
observing epoch, a window of span up to 15 days is de- 
fined over which an independent DM is fit for. Epochs 
for which no dual- frequency data exist within a 15 day 
range were excluded from the analysis. This leads to a 
DM versus time measurement for each pulsar, as shown 
in Figure O As previously noted, PSRs J1853-hl308, 
J1910-I-1256 and B1953-h29 have mostly single-band ob- 
servations, so the DM variation has not been modelled 
for these three sources. 

During the course of this analysis, we discovered addi- 
tional radio frequency-dependent trends in pulse arrival 
times that were not well described by the ly^^ dispersion 
relation. We attribute these to the intrins ic evolution of 
the p ulse profile shape with frequency fe.g. lKramer et all 
Il999| ). The interstellar medium is another possible cause, 
however since the effect does not appear to be time- 
variable an ISM explanation seems less likely. Regard- 
less of its physical origin, any profile shape change versus 
frequency can lead to systematic TO A biases as follows: 
Although we have used separate template profiles per re- 
ceiver, within each receiver band the template profile is 



constant. If the true profile shape is changing versus fre- 
quency, this will cause small but measureable systematic 
effects. To correct for these biases, we have included as 
free parameters in the timing fit a constant (in time) off- 
set - also known as a "jump" - for each 4 MHz frequency 
channel. 

A summary of the timing analysis and results is pre- 
sented in Table [2l For each pulsar, the total number of 
TOAs, and the total number of fit parameters is given. 
These are divided into those relating to the DM vari- 
ation, frequency-dependent terms, and the other stan- 
dard spin, astrometric and binary parameters. The suc- 
cess of the timing model fit can be characterized by an- 
alyzing the post-fit residuals, i.e. the difference between 
the observed and model-predicted arrival times. The 
uncertainty- weighted root-mean-square (RMS) residual 
value and normalized values as determined directly 
from the fit are listed in the table. The majority of 
values fall near 1, but are as high as ~4 in some cases. 
These values have been computed without the use of any 
multiplicative or additive modifications to the TOA un- 
certaintiesI3 Traditionally pulsar timing results have 
been presented using TOAs which come from an aver- 
age of all data taken during a given day. Here, our fit 
procedure requires the multifrequency data be kept sep- 
arate, so we cannot perform this averaging pre-fit. To 
facilitate comparison with previous work, the final three 
columns in Table [5] show the "epoch-averaged RMS" 
computed by first averaging the raw post-fit residuals 
down to one point per receiver per epoch, then taking 
their weighted RMS. The RMS from the low- frequency 
and high-frequency bands for each pulsar are shown sep- 
arately as well as the combined value using all the data. 
Due to the DM{t) fit, the low-frequency RMS tends to be 
suppressed, making it less useful as a simple characteriza- 
tion of the timing. The two best pulsars in the set, PSRs 
J1713-H0747 and J1909-3744, both have epoch-averaged 
RMS in the ~30-50 ns range with comparable values for 
both the high-frequency and combined versions of this 
statistic. 

4. GRAVITATIONAL WAVE ANALYSIS 

The presence of gravitational waves (GW) along the 
line-of-sight from a pulsar to Earth alters the effective 
path length in a time-variable manner, result ing in ex- 
tra p erturbations in pulsar timing residuals (jDetweiled 
|1979() . The observed amplitude of a single pulsar's tim- 
ing residuals can therefore be used to limit the strength 
of an y GW that may exist (jKaspi et al.lll994l : iJenet et all 
120061 ). However, as an observed timing perturbation 
could also arise from a number of non-GW sources, this 
method can not be used to definitivel y confirm the pres- 
ence of GW. It was first noted by iHellings fc Downi 
(1983) that a GW signal induces correlated variations 
in the timing residuals of a set of pulsars. The form of 
this correlation is unique to GW among the expected 
perturbations, and attempting to detect it is the ba- 
sis for most current pulsa, r timing array (PTA ) efforts 
(jvan Haasteren et al.ir201lHYardley et al.ll2011| ). includ- 
ing this work. 

24 I.e., the "EFAC" and "EQUAD" options respectively in 
TEMPO. These parameters were sometimes used in previous tim- 
ing analyses to compensate for unexplained systematic errors in 
TOAs or their uncertainties. 
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Fig. 3. — Timing summary for PSR J1713+0747. The top panel shows residuals from the multi-frequency TOAs used in the timing fit. 
The middle panel shows the same residuals averaged down to one point per band per day. The bottom panel shows the measured variation 
in DM as a function of time. 



In this analysis, we will search for a stochastic grav- 
itational wave background (GWB) signal in the pulsar 
timing results presented in We assume the GWB will 
take the form which has become standard in this field 
- a power-law frequency spectrum and isotropic angu- 
lar distribution. This signal is expected to be generated 
by the sum of unresolved supermassive black hole binary 
systems with masses of ^ 10* Mq and orbital periods of 
1-10 years. In this case, the characteristic strain spec- 
trum is expected to have a "red" power-law spectral in- 
dex a = -2/3 fe.g.. iJaffe fc Backed [200l iSesana et al.l 
120081) ■ with a possible br eak near 10 nHz du e to the fi- 
nite number of sources (ISesana et aUMM ). A GWB 
of this form could also be ge nerated from cosmic su- 
perstrings, with a — —7/6 (D amour fc VilenkinI 120051 : 
[Siemens et al.ll200'^. or as infiation arv "relics" with spec- 
tral index a = -1 ()Grishchukll2005[) . Any GW signal will 
produces correlation in the timing fluctuations of pairs 
of pulsars. In the specific case of an isotropic GWB, 
the amount of correlated power is a function only of the 
angular separation of the two pulsars in the pair, and 
has a characterist i c fun ctional form first predicted by 
IHellines fc Downg (flQSl . 

In this paper, we adopt definitions of the expected 
gravitational wave spectrum and its effect on tim- 
ing consistent with previous papers on the topic (e.g., 
IJenet et al.l I2006t Ivan Haasteren et alj 120111 ) . In partic- 
ular, we assume a power-law spectrum in characteristic 
strain. 



hcif) = Af, f 



(1) 



Here, Aj^ is the unknown GW spectrum amplitude at 
a reference frequency /q. For consistency with previous 
literature, we set /o = 1 yr~^, and will call the resulting 
amplitude Ai. This GWB produces a fluctuation y{t) in 
the pulse times of arrival from a given pulsar, with power 
spectrum given by 



Syif) 



1 



127r2 f 



hcif)' 



(2) 



It is important to note that y{t) represents the pre- fit 
contribution of the GW signal to the pulse TOAs. The 
effect of the timing model fit will be considered in §4.11 
Also important is that t his formulation of S ,, (f) is consis- 
tent w ith t hat used by IJenet et al. (2006) , H obbs et al.l 
(j2009() , and Ivan Haasteren et all ( 2011 but is a factor 
of 3 smaller than that used bv IJenet "etan (|2005h . This 
results in a factor of \/3 difference in limits on Ai de- 
pending on which definition is in use, and care should 
therefore be taken when directly comparing the various 
pubhshed limits. 

It is useful to compute the expected time-domain cor- 
relation between pairs of timing fluctuations from pulsars 
a and b. 



c: 



(ab) 



E {ya{U)ytitj)} = CyiU - tj)a0ab). (3) 



-^ ^ While the e quations given bv lvan Haasteren et al.l I I2011I1 used 
the IJenet et al.l |(2005J) definition, their published limit of Ai < 
6 X 10~^^ was computed using the same scaling as our Eqn.[2](van 
Haasteren 2012, private communication). 
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Here, Oat is the angular separation between pulsars a and 
b, and C(^a6) is the Hellings-Downs function describing 
the expected angular correlation for a isotropic GWB. 
Cy{T) is the GW signal autocorrelation (Fourier trans- 
form of Sy{f)), and can be computed analyticall y for a 
power - law sp ectrum as presented bv .van Haasteren et alJ 
i|2009l 120111 ). Additionally, here and in the following, 
E{-} represents the statistical expectation value, or av- 
erage over many realizations of the enclosed quantity. 
This description assumes the GW signal foll ows wide- 
sense stationary statistics (e.g.. iPapoulis fc^ illai 2002). 

In the following sections we first present methodology 
for analyzing the timing residuals from 331 a-nd discuss 
the effect that fitting the timing model has on the statis- 
tics of post-fit residuals. We will then apply these meth- 
ods to first determine an upper limit to the stochastic 
GWB amplitude using a single pulsar, and then attempt 
to detect or limit the GWB by looking for angular cor- 
relations. 

4.1. Effect of the Timing Model Fit 

In order to use pulsar timing data to detect gravita- 
tional waves, one must account for the fact that the pul- 
sar parameters (spin period, astrometric parameters, etc) 
are not known a priori, and need to be determined from 
the same data that are used for GW detection. Especially 
in the case of red GW spectra, a large fraction of the GW 
signal power is covariant with the "long-term" pulsar pa- 
rameters such as spin period and spin-down rate, and 
cannot be unambiguously separated from these intrinsic 
pulsar features. Previous analyses have typically dealt 
with this by decomposing timing residuals usin g a set 
of po l ynomials of cubic and higher order ( Kaspi et al.l 
et al.l |2006() . These are approximately or- 
thogonal to the rest of the timing fit, and capture most 
of the low-frequency GW powe r. Recent analyses by 
Ivan Haasteren et all ()2009l 1201 ID use a Bayesian frame- 
work to limit or detect GW signals in the d ata, marginal- 
izing over the timing model parameters. lYardlev et al.l 
(|20ll use an extensive series of simulations to charac- 
terize the effect of timing parameter fitting on the GW 
signal in the frequency domain. 

In this work, we perform the GW analysis in a sep- 
arate step from the timing fit, while using the proper- 
ties of the timing fit to both determine how much GW 
power is absorbed by it, and how to best detect or limit 
that which remains. This is built around a computa- 
tion of the statistics (covariance matrix) of the post-fit 
res iduals. Th is anal ysis draws heavily on that presented 
by iDemores^ (|2007( ). which can be referred to for ad- 
ditional discussion. Further detailed discussi on and ap- 
plicat ion of this method will be presented by lEllis et all 
(|2012f ). It is also interesting to note that a very similar 
calculation app ears in the B ayesian method presented by 
[van Haasteren et all (120091 ) . as a means of marginalizing 
over the timing fit "nuisance" parameters. 

To begin we note that the fit that determines the pul- 
sar parameters is a weighted least-squares minimization, 
which in matrix terms is the solution of the normal equa- 
tions (see for example iPress et al]|1992[ ): 

A^WAa = A^Wy (4) 
Here y is the vector of input data (sometimes called "pre- 



fit residuals" ) , A is the fit design matrix whose columns 
are the fit basis functions, a is the vector of fit parame- 
ter values, and W is the weighting matrix. In a standard 
weighted least squares fit such as is used in this paper, 
W is a diagonal matrix whose entries are calculated from 
the TOA uncertainties, Wa = cr~^. However, we write 
Eqn. [4] in this general form to emphasize that the fol- 
lowing methods are equally valid for a generalized least- 
squares fit in which W can have non-zero off-diagonal 
elements. The generalized approach has recently been 
proposed for pulsar timing analyse s as a way of hand ling 
correlated noise in the timing data ()Coles et al.ll201l|) . In 
the generalized case, W is the inverse of a non-diagonal 
noise covariance matrix. 

The post- fit residuals r = y — Aa = Ry are therefore a 
linear function of the input data, and can be determined 
by applying the residual projection operator R: 

R = I - A (A^WA) A^ W (5) 

A key feature of writing the analysis in this form is that 
althought R does depend of the weighting matrix W, it 
does not depend on the data values y. Therefore, for 
a given pattern of weights (or uncertainties) , the action 
of the fit can be studied independent of a specific data 
realization. This feature can be used to eliminate depen- 
dence on simulation in interpreting results and determin- 
ing GW limits, and can provide additional insight into 
the properties of the fit procedure. In particular, given 
a known or assumed input data covariance matrix Cy, 
the resulting residual covariance matrix Cr can be easily 
calculated: 

Cr^E {rr'^} = HE {yy^} R^ = RCyR^ (6) 

The expected cross-covariance between the timing resid- 
uals of a pair of pulsars (a and b) can be calculated sim- 
ilarly: 

Cr^,rt, = E {rarb^} = RaCy^,yj,Rb^ (7) 

It is worth noting that the ensemble of post-fit residuals 
is in general non-stationary. That is, the elements of the 
covariance matrix, Cr,ij , are not simply a function of the 
time lag = tj —ti, even if the pre- fit stochastic process 
represented by Cy does have this property. Therefore the 
action of the timing fit cannot be fully represented as a 
filter or frequency domain transfer funct ion, as has some- 
times been proposed in the past (e.g., IBlandford et al.l 
1984; Jcnet et al. 2005). This feature of the statistics of 
the residuals, along with the irregular sampling, varying 
data span and quality, and presence of steep-spectrum 
red noise processes, present significant obstacles to a fre- 
quency domain analysis of real- world pulsar timing data. 
This motivates our decision to perform the GW analysis 
in the time domain. 

4.2. Single-pulsar Analysis 

We can use the methods discussed in the previous sec- 
tion, along with timing data from a single pulsar, to com- 
pute upper limits on the strength of the GWB. As previ- 
ously discussed, there is no definitive way to distinguish 
GW-induced timing fluctuations observed in a single pul- 
sar from those due to intrinsic pulsar spin irregularity or 
other non-GW effects. However, in the low-amplitude 
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GW regime, focusing on the single best pulsar in this 
manner can provide some of the most constraining up- 
per limits. 

Given a GW spectrum of assumed power- law shape but 
unknown amplitude, we compute the pre-fit GW covari- 
ance matrix up to an ove rall scaling by Af (see 

Eqn. El and Ivan Haasteren et all (|2009f )). This is then 
converted to the residual covariance matrix using 
Eqn. \6\ Diagonalizing provides an orthonormal 

basis of eigenvectors that can be used to decompose the 
timing residuals. This basis represents the GW sign al us- 
ing the smallest possible number of components and is 
completely orthogonal to the timing fit. In other words, 
this basis optimally captures the portion of the GW sig- 
nal in the data that is not absorbed by the timing model 
fit. Transforming the residuals into this basis produces a 
set of projection coefficients c^. The associated eigenval- 
ues (Ai) give the relative level of GW signal power in each 
component, and by assuming a nominal value of Ai, can 
be scaled to produce expected mean-square timing resid- 
ual values. Figure |4] shows the measured and expected 
component amplitudes for the two best pulsars in our set, 
PSRs J1713-H0747 and J1909-3744, assuming a = -2/3 
and Al = 10~^^. In these plots and the following discus- 
sion, the residual coefficients and eigenvalues have been 
normalized to represent the RMS residual due to each 
component. The quadrature sum of all coefficients for a 

1 /2 

given pulsar c^) reproduces its full RMS residual 
as presented in Table [2] 

It is immediately apparent that the data shown in a 
plot such as Figure |4] can be used to place a limit on the 
strength of a power-law GWB. Roughly speaking, any 
GW signal present in the data must lie near or below the 
observed data values. For measured component values Ci 
and eigenvalues Ai, we combine several components into 
a single limit by forming the likelihood function for the 
parameters Ai and cr, the level of white noise present in 
the data: 

logi(Ai,a) = - i E {log27r(A?A, + a') + ^2^^] 

(8) 

This formulation assumes both the GW signal and the 
white noise follow Gaussian statistics; in this case each 
component value Ci is drawn from a normal distribution 
with zero mean and variance equal to A\Xi+u^ . The dif- 
ferent components are by construction independent due 
to the choice of basis. The product of all the compo- 
nent distributions for a given pulsar results in the like- 
lihood function of Eqn. [S] Using uniform/positive pri- 
ors on a and Ai, and marginalizing over cr, we convert 
the likelihood function above into a posterior probabil- 
ity distribution for Ai, shown in Figure [5] For the data 
from J1713-f0747, 95% of the distribution faUs under 
Al < 1.1 X 10~^^ for an a = —2/3 spectrum. This is the 
most constraining single pulsar in our set. Similar upper 
limits computed for each of the pulsars at a variety of 
astrophysically relevant a are listed in Table |31 

^® This is very similar to prinicpal components analysis (PC A). 
The only difference is that PCA is typically based on an empirical 
data covariance matrix rather than an assumed covariance matrix 
as we have used here. 



While we have framed the discussion so far in terms 
of a GWB signal, this analysis method is more funda- 
mentally a way to quantify timing perturbations of a 
known/assumed spectral shape that may exist in the 
data, regardless of their cause. In this context we can 
look at these single-pulsar results as a test for "red" 
power-law timing noise in the data set. In addition 
to the 95% upper limit on Ai already described, two 
other statistics are useful: The maximum-likelihood es- 
timate Al and the ratio R = L{0,(Jo)/L{Ai,a) between 
the maximum likelihood values obtained by either fix- 
ing ^1 = or allowing it to vary. Ai characterizes the 
strength of the red noise, and R characterizes its signifi- 
cance in the data. Values of R near 1 indicate consistency 
with a white-noise-only model, while very small values of 
R indicate significant red noise is detected. These three 
values are listed for all pulsars for a selection of vari- 
ous expected GW spectral indices in Table [31 We see 
that two of the pulsars - J1643-1224 and J1910-hl256 
- show very significant red noise (i? < 0.01), and two 
others - J1640-1-2224 and B1953-h29 - show mildly sig- 
nificant red noise (0.01 < i? < 0.1). The remaining 13 
sources are consistent with white noise only. Additional 
detailed analyses of the noise in these data are ongoing, 
and include both a Bayesian ana lysis covering a wider 
range of power-law s pectral index (lEllis et al.ll2012l ) , and 
an application of the 'Coles et al.' (j2011|) methodologv to 
our data set (Perrodin et al. 201^. 

We will not attempt to definitively determine the cause 
of the red noise in the four sources noted here. However 
we speculate that at least for three of them a likely source 
is the interstellar medium, rather than processes intrinsic 
to the pulsars. As explained above, J1910-I-1256 and 
B1953-I-29 do not include DM(t) fits so their timing wiU 
be affected at some level by unmodelled DM variation. 
J1643— 1224 has the second highest DM of any pulsar in 
our set (after B1953-I-29), and has a high predicted level 
of interstellar scatter-broadening, suggesting multipath 
ISM effects may be an issue here. 

4.3. Verification 

The methods discussed in the previous section do not 
rely at all on simulated data to produce a GW upper 
limit or red noise estimate. However, it is still a valuable 
exercise to analyze a simulated GW signal using this ap- 
proach, as a test of its validity. To do so, we added 
simulated GW to our timing data for J1713-I-0747 , using 
the TEMP02 GWbkgrd plugin (Hobbs et al.l 12009V 
GW were added using a = —2/3 and Ai = 10""'^'*. These 
data were then analyzed using the eigenvalue procedure 
described above. Figure [5] shows the root-mean-square 
component values (q) from 1000 GW realizations, along 

1/2 

with the square root of the associated eigenvalues (A, ' ) 
scaled to the known GW signal level. The excellent 
agreement between these two lines shows that our al- 
gorithm is correctly predicting the amount of GW power 
present in the timing residuals. Figure [7] presents results 
of computing maximum-likelihood and 95% upper limit 
amplitude values from the simulated data. These statis- 
tics agree well with the known amount of GW signal in 
the simulated data (grey region in plots). 

As an additional test of this method, we analyzed the 
pubhcly-available seven-year B1855-f09 timing data set 
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Fig. 4. — Eigenvalue spectra for J1713+0747 (top) and J1909— 3744 (bottom). In each plot, the line labelled "Data" shows the absolute 
value of the residual coefficients Ci when transformed into the optimal basis for an a = —2/3 power-law GW spectrum. The "Exp. GW" 

line shows the expected amplitude (square root of the eigenvalues, aV'^) for an Ai = 10~^^ GW signal level. The fact that the data lines 
are flat reflects the fact that the residuals are dominated by white noise. 



TABLE 3 

Results from single-pulsar GW analyses. 



Source 


logio R 


a = -2/3 
Max-L Ai 95% 
(10-15) 


Ai limit 
(10-15) 


logio R 


a = -1 
Max-L Ai 

(10-15) 


95% Ai limit 
(10-15) 


logio R 


a = -7/6 
Max-L Ai 95% Ai limit 
(10-15) (10-15) 


J0030-I-0451 


0.00 


0.0 


76.5 


0.00 


0.0 


47.8 


0.00 


0.0 


37.8 


J0613-0200 


0.00 


0.0 


100.9 


0.00 


0.0 


63.1 


0.00 


0.0 


49.9 


J1012-I-5307 


0.00 


0.0 


153.8 


0.00 


0.0 


109.7 


0.00 


0.0 


93.5 


J1455-3330 


0.00 


0.0 


242.8 


0.00 


0.0 


162.6 


0.00 


0.0 


133.1 


J1600-3053 


-0.00 


58.9 


973.9 


-0.01 


62.6 


1173.6 


-0.03 


68.1 


1240.4 


.11640+2224 


-1.31 


68.1 


501.6 


-1.84 


40.3 


413.3 


-1.90 


28.1 


350.2 


J1643-1224 


-4.57 


182.9 


455.3 


-4.57 


122.5 


333.6 


-4.70 


99.5 


284.6 


J1713+0747 


-0.22 


1.9 


11.3 


-0.22 


1.1 


7.6 


-0.29 


0.8 


6.2 


J1744-1134 


0.00 


0.0 


184.1 


0.00 


0.0 


119.1 


0.00 


0.0 


96.2 


J1853+1308 


-0.45 


23.6 


131.2 


-0.36 


8.3 


86.7 


-0.35 


5.3 


70.0 


B1855+09 


-0.31 


14.0 


70.0 


-0.59 


8.6 


45.6 


-0.65 


6.4 


36.3 


J1909-3744 


0.00 


0.0 


39.4 


0.00 


0.0 


26.8 


-0.00 


1.0 


25.0 


J1910+1256 


-4.66 


41.7 


229.7 


-4.91 


23.6 


144.6 


-4.85 


17.4 


115.1 


J1918-0642 


-0.31 


64.8 


545.0 


-0.21 


38.9 


383.1 


-0.11 


10.5 


282.6 


B1953+29 


-2.11 


274.9 


1877.7 


-1.43 


188.0 


1249.0 


-1.52 


160.4 


1110.5 


J2 145-0750 


0.00 


0.0 


568.0 


0.00 


0.0 


372.6 


0.00 


0.0 


296.6 


J2317+1439 


-0.20 


43.4 


383.1 


-0.13 


26.4 


195.9 


-0.07 


20.9 


148.6 


B1855+09 (KTR94) 
B 1855+09 (combined) 


-0.27 
-8.53 


5.6 
5.0 


31.3 
13.1 


-0.15 
-7.09 


1.8 
1.2 


19.3 
6.5 


-0.12 
-7.41 


1.1 
0.7 


15.2 
4.8 
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GWB Amplitude (A^) 



Fig. 5. — Probability distribution for the GWB amplitude Ai based on J1713+0747 timing, and assuming ct = —2/3. 95% of the 
distribution is contained in Ai < 1.1 X 10~^^ (vertical line), the maximum-likelihood value A\ = 1.9 X 10~^^, and the likelihood ratio 
between this point and Ai = is i? = 0.6. In this case, a non-zero value of A\ provides the best fit to the data, but without much statstical 
significance over a white noise only model. These statistics for all sources, considering several different values of a, are presented in Table [3l 
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published bv lKaspi etall ()1994 hereafter KTR94), from 
observations taken during 1986-1992. This data set has 
been used for several previous GW analyses, so provides 
a useful benchmark for testing different analysis meth- 
ods. For this comparison we focus on the GW spectrum 
originally considered by KTR94, with a = —I. From the 
KTR94 data set alone we find no significant evidence for 
red timing noise (see results in Table [3]) , and using our 
method it sets a 95% upper hmit of Ai < 1.9 x 10^^^. 
This is somewhat more conservative than the published 
limits of ^1 < 1.3 X 10-1'^ (KTR94) and Ai < 1.4x lO^i^ 
(jJenet et al.ll2Q06l ) obtained from the same dataEH This 
difference may be due to the varying strategies used by 
these three analyses to account for effects of the timing 
model fit and the possibility of existing red noise in the 
data. 

Combining the KTR94 data with our newer observa- 
tions of this pulsar into a single analysis produces a limit 
Ai < 6.5 X 10~^^, comparable to what we get from our 
J1713-f0747 data alone. However, very significant red 
noise is detected in the combined B1855-I-09 dataset (Ta- 
ble [31) . A fundamental feature of red noise processes is 
that their amplitude grows with increasing data span. 
Therefore this detection in the combined data set is not 
inconsistent with the lack of such a detection in either of 
the two data subsets on their own. Another contributing 
factor may be unmodeled variation in DM. The KTR94 
data were obtained at a single frequency and as such can- 
not be corrected for DM variation. In our newer data, a 
significant DM(t) trend is clearly visible for this pulsar 
(Figure m]). 



4.4. Cross- correlation Analysis 

To detect a GW signal present in the timing data, we 
must consider correlations among the timing residuals of 
different pulsars, as previously discussed. This is compli- 
cated by the fact that residuals from each pulsar are pro- 
duced by separate timing model fits. The fits potentially 
incorporate a wide variety of effects (e.g. binary param- 
eters) that vary from pulsar to pulsar. Furthermore the 
span and cadence of observations is not necessarily iden- 
tical between all the sources, and each may be affected 
by a different level of red timing noise from either GW 
or other sources. Correctly measuring the correlation 
between two such residual timeseries is therefore more 
complicated than simply multiplying and averaging the 
two signals, as would be done in a standard correlation 
analysis. 

Referring to the discussion in t ^4.1l and in particular 
Eqn. [71 given the actual timing model fits used for any 
two pulsars and an expected gravitational wave spec- 
trum, we determine the expected cross-correlation for 
the GW signal between all pairs of residuals of each pul- 
sar pair (C^^ ). We use the elements of this matrix. 



along with the estimated noise parameters from t ^4.2l as 
weights for a cross-correlation sum to estimate the GW 



Limits originally published in units of ^g-uih? were converted 
to equivalent Ai using Eqn. 3 of lJenet et al.l 120061) . 



signal power correlated between pulsars a and b: 

l^zjkl^z C-- j.„ Sfc \^ )kl 



Pab 



The components of this expression are: 



(9) 



• The sum indices i and j run from 1 to Na, the 
number of TO As for pulsar a. Similarly, k and I 
run from 1 to Nh. 

• rj-"-* and r^-* are simply the post-fit timing residuals 
for pulsars a and b respectively. 

• C'*°**-°-' and C'^i*'^''^ are the total post-fit covariance 
matrices for each pulsar. These are determined us- 
ing the maximum-likelihood red noise parameters 
found in the earlier single-pulsar analysis, com- 
bined with the TOA measurement uncertainties, 
and modified by the timing model fit (see Eqn. ^ : 

Cj°* = R(ii'cGW + W-i)R^ (10) 



C, 



GW{a,b) 



J). are the elements of C^^, the expected 

post-fit cross-covariance matrix for the GW signal 
in pulsars a and 6, computed using Eqn. [3 and the 
assumed power-law spectral shape. 

The presence of the inverse covariance matrices in Eqn. [SI 
effectively serves to "whiten" the timing residuals before 
computing the cross-correlation sum. This reduces the 
weight of data from pulsars that have red noise compo- 
nents, potentially including any "self- noise" contribution 
from the stochastic GWB itself. This formulation is con- 
ceptually very simila r to the freque n cy-do main optimal 
statistic presented bv lAnholm et al.| (|2009|) . 

The post-fit covariance matrices C*°* are singular, due 
to the degrees of freedom removed by the timing fit, or 
equivalently by the application of the R operator. There- 
fore the matrix inverse technically does not exist. In 
this situation it is appropriate to use a singular value 
decomposition-based pseudoinverse in its place. This is 
equivalent to performing the correlation in the subspace 
of signals that are orthogonal to (i.e., not absorbed by) 
the timing fit. The uncertainty on pab is given by: 



^GW{a.b) /^tot(b) 



a 



GW{a,b) 



kl 



\ijkl 

The resulting pab for all 136 pulsar pairs, with a — 
—2/3, is shown in Fig. [8l We search for the presence of 
the GWB in these data by fitting them to an amplitude 
(proportional to Af) times the Hellings-Downs angular 
function ({Oab)- The best fit, also shown in Fig. HI results 
in Aj = (-10 ±26) X 10-30 for a = -2/3. This is consis- 
tent with no detectable GWB in the data, and also can 
be interpreted as a 2-cr upper hmit of Ai < 7.2 x 10"^^, 
an improvement over the single-pulsar limits presented in 
^A.2^ The reduced-x^ of this fit is 0.95, providing addi- 
tional confidence that the uncertainty estimate of Egn. lTll 
is correct. 
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Fig. 6. — Eigenvalue spectra for J1713+0747 with simulated GW injected at an amplitude of Ai = 10"^''. As in Figure|4l the "Exp. 
GW" line shows the expected GW amplitude determined from the eigenvalue analysis, scaled to the known signal level. The "Data" line 
shows the root-mean-square residual in each component, averaged over 1000 realizations of the simulation. The correct GW signal level is 
observed in the data, until the level drops below the white noise near component 7. 



A table of the top 15 correlation measureraents, sorted 
by increasing uncertainty, is presented in Table |4l It is 
clear that the measurement is dominated by pairs in- 
volving one or both of the two best-timing pulsars in the 
set, J1713-K0747 and J1909-3744. In total, 8 separate 
pulsars contribute to these 15 lowest-uncertainty points. 
The limit being set primarily by a small number of pul- 
sars is due to the fact that our current data set is in a 
white noise dominated, low GW amplitude regime. As 
the length and quality of such datasets increase, they 
are expected to become increasingly dominated by red 
noise components, either from the GWB itself or from 
intrinsic pulsar timing noise. In the red noise dominated 
and/or strong GW regime, a large number of pulsars 
(likely ~20-50) is required to improve the significance of 
a GW detection penet et al.l l2005t iCordes fc ShannonI 
1201 It) . This provides motivation for continuing to ob- 
serve many MSPs, even if upper limits are currently only 
dominated by a small subset of the pulsars. 

Table m also lists the best-fit and its uncertainty for 
three different spectral indices. The resulting 2-a upper 
limits are Ai = 7.2 x lO-^^^ 4.1 x 10'^^, and 3.0 x lO^^^ 
for a — —2/3, —1 and —7/6 respectively. The change 
in value of Ai with a is primarily due to the fact that 
although we are using 1 year~^ as the reference GW fre- 
quency, the measurement is most sensitive to frequencies 
near T~^, where T ~ 5 years is the effective length of the 
multi-pulsar data set. Therefore, we expect the various 



limits on Ai to scale with a as follows: 

A,{a)=AT(-^y (12) 
VI year/ 

Using our three measured Ai limits to empirically de- 
termine At and T gives At = 2.26 x 10~" and T = 
5.54 years (see Figure [5]). These values can be used 
togther with Eqn. [T^] to convert our measurement into 
an approximate limit on Ai for values of a that we have 
not considered here. 

5. DISCUSSION 
5.1. Astrophysical Consequences 

As mentioned previously, the strongest expected source 
of GW in the nanohertz band is the ensemble background 
from supermassive black hole binaries. It is on these 
which we will focus our discussion here. The GWB pre- 
dicted from the supermassive black hole binary popu- 
lation is dependent on a number of observationally ill- 
constrained values: e. g. the ubiquity of supermassive 
black holes in galaxy centers, the prevalence of merging, 
as well as the timescale and dominant mode of black hole 
growth (via ac cretion or merge r with other black holes). 

Recently, Se sana et al.l (|2008l ) considered a broad range 
of probable parameter uncertainty values applied to four 
leading scenarios for black hole growth in the context 
of hierarchical merging. Of most importance here, they 
found that when accounting for uncertainties in all model 
parameters, the estimated range of predicted Ai values 
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Fig. 7. — Histograms of maximum-likelihood GW amplitude (upper panel), and 95% upper limit (lower panel) computed from 1000 
realizations of a simulated Ai = 10~^^, a = —2/3 GWB added to J1713+0747 timing data. The shaded range repres ents the true GW 
signal level, including both the simulated signal and any existing GW signal in the data, up to the limit presented in t|4.2l The vertical 
line in the lower panel shows the 5% quantile of the histogram values. 



TABLE 4 

Cross-correlated power measurements and GW results. 



Pulsar a 


Pulsar b 


Angle Oab 
(deg) 


a = -2/3 
Cross-power Uncertainty 
Pat (10-30) a,„, (10-30) 


a = -1 
Cross-power Uncertainty 

Pat (10-30) ap„, (10-30) 


a = 

Cross-power 

Pat (10-30) 


-7/6 

Uncertainty 
-P„. (10-30) 


J1713-I-0747 


J1909-3744 


53.0 


-28.4 


13.2 


-8.8 


4.4 


-5.9 


2.7 


J1713-I-0747 


J1744-1134 


20.8 


-13.1 


23.4 


-3.5 


7.7 


-1.6 


4.0 


J0613-0200 


J1713+0747 


164.0 


24.0 


33.4 


7.1 


11.0 


3.8 


5.9 


J1012-I-5307 


J1713+0747 


92.8 


4.8 


35.0 


1.1 


11.2 


0.6 


6.0 


J1744-1134 


J1909-3744 


32.4 


8.3 


40.7 


-0.6 


12.6 


0.4 


8.0 


J0030-I-0451 


J1713+0747 


108.2 


21.2 


41.0 


6.7 


12.6 


3.8 


6.5 


J1713-I-0747 


B1855+09 


25.7 


-44.0 


49.5 


-11.5 


17.9 


-5.5 


9.9 


J0613-0200 


J1909-3744 


138.2 


-40.9 


59.4 


-13.0 


18.3 


-7.8 


11.7 


J1012-I-5307 


J1909-3744 


145.2 


-9.9 


62.6 


-4.8 


18.7 


-3.1 


11.9 


J0030-I-0451 


J1909-3744 


85.3 


-44.4 


75.0 


-14.4 


21.2 


-8.3 


13.0 


J1713-I-0747 


J1853+1308 


25.2 


-71.8 


82.1 


-18.7 


19.5 


-9.8 


9.6 


B1855+09 


J1909-3744 


47.5 


110.2 


97.3 


37.1 


33.1 


20.0 


21.2 


J0613-0200 


J1744-1134 


164.6 


31.6 


103.2 


9.3 


31.3 


5.2 


16.8 


J1012-I-5307 


J1744-1134 


113.0 


77.4 


108.6 


21.7 


32.2 


11.5 


17.1 


J0030-I-0451 


J1744-1134 


102.2 


39.9 


140.0 


11.6 


38.2 


6.4 


19.4 


Best-fit A'j (10"^" 


) 


-10 ±26 


-3.7 ±8.4 


-1.9 ±4.6 
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Fig. 8. — Measured cross-correlated power p^i, as a function of separation angle O^b for pairs of pulsars in our set, with error bars showing 
1-<T uncertainty. Power is normalized relative to an Af = 10-3°, a = -2/3 GWB. The lines show the ±2 cr fit to the amplitude of the 
Hellings-Downs function CW- All 136 cross-correlation points were used for the fit, however for clarity the 15 lowest-uncertainty values are 
denoted with solid/bold symbols. 
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Spectral index (a) 

Fig. 9. — Measured 2-a upper limits on A\ as a function of GW 
spectral index a (squares). The values are consistent with a simple 
scaling based on At = 2.26 X 10"^* at T = 5.54 years and Ean.[T2l 
(line). This relationship can be used to convert these results to 
equivalent limits for other values of a. 

produced a relatively narrow predicted amplitude range 
of 1 X 10-16 < ^1 < 3 X 10^1^ While we are not yet 
sensitive enough to impact even the upper boundary of 
their predicted range, it is pertinent to note that our 
current sensitivity is approximately a factor of 2 away 
from the expected highest prediction. If adequate GW 
sensitivity can be achieved, pulsar timing will probe as- 
trophysical quantities not readily accessible by electro- 
m agnetic observ a tions. The most uncertain factors noted 
bv lSesana efall (|2008D were the supermassive back hole 
mass function, the local merger rate of massive galax- 
ies, the amount of merger-induced accretion onto black 
holes, and the post-merger black hole inspiral rate. A 
measurement or constraint on the mass function leads 
to limits on galactic host/black hole mass relationships, 
particularly in the most massiv e regime. 

The other relevant finding of ' Sesana et al.l (|2008D was 
that with a paucity depending on the input model, GW 
spectral bins above ~1 — 5 x lO"* Hz were populated by 
increasingly few black holes - at the highest frequencies, 
sometimes one or none. A lack of black holes contribut- 
ing to the high-frequency GW spectrum could affect us 
in two ways. It c ould change the ob served GWB spec- 
tral slope. In fact, Sesana et al.l ([2008) uses an expanded 
version of our Eq. 1 to include a break in the spectrum. 
Alternately it could induce non-stochastic signals, requir- 
ing a directional/point-source analysis approach to de- 
tect the emitted signal. Either of these consequences 
would corrupt the sensitivity and applicability of the 
analysis performed here, particularly if the analysis were 
performed on data spans of ^3 years. All but one of the 
pulsar data sets presented here exceed this length, thus 
the influence of this issue is expected to be minimal. As 
the contributing population is expected to increase at 
lower frequencies, the impact of this problem will de- 
crease as further data are collected on these pulsars, and 
as we reach the aforementioned sensitivity improvements 
that come with increased data span. Nonetheless, further 
studies must be performed to address the effects of the 
A^-source contribution issue, particularly to understand 
its effect on the expected Hellings-Downs function. 

5.2. Cosmic Strings 



Cosmic strings or superstrings are linear-dimensional 
structures proposed to arise due to phase transitions 
in the early universe. Interactions between mem- 
bers of a network of such objects may create cusps 
or loops that are precicted to radiate significant GW, 
producing a stochastic GWB wi th a ~ —7/6 in 
the nanohertz frequency band (e.g., iDamour fc VilenkinI 
120051: [Siemens et al.ll2007[) . The amount of GW produced 
depends on physical parameters such as the string loop 
size scale, the reconnection probability and the string 
tension. Conversely, measured GW limits constrain the 
possi ble values of these pa rameters. Following the analy- 
sis of lSiemens et al.l (|2007f) , our results for cosmic strings 
take the form of sections of cosmic string parameter space 
that are ruled out or allowed by our upper limit on the 
gravitational wave background. 

In the case where cosmic string loop sizes are set by 
gravitational back reaction ("small loop" case) we ex- 
plore the three-dimensional parameter space of cosmic 
string dimensionless tension G/i, reconnection probabil- 
ity jp, and loop size s parametrized by e (for details, see 
ISiemens et al.ll2007D . The shaded areas of Fig. [TU] show 
the regions of the e-Gfj, plane ruled out for four values 
of the reconnection probability p = 1, 10^^, 10"^ , and 
10~^. For example, for e = 1, and p — 10^^ the string 
tension G/x is constrained to be less than about 2 x 10"^^. 

If the size of cosmic string loops is instead given by the 
large scale dynamics of the n etwork, as suggested by re - 
cent numerical simulations ( Blanco- Pillado et al.ll2011|) . 
we fix the cosmic string loop length and explore the two- 
dimensional parameter space of reconnection p r obabi l- 
ity and string tension following ISiemens et al.l ()2007D . 
These results are shown in Fig [TTJ The stochastic back- 
ground produced in these models is substantially larger 
therefore our constraints on parameters are correspond- 
ingly tighter. For example for a reconnection probability 
p = 1, Gfi < 10"^, and for a reconnection probability 
of p = 10^"^ all cosmic string tensions above 10^"'^^ are 
ruled out. 
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Fig. 10. — Cosmic string parameter space constraints from our 
measurement, in the small loop case. The shaded areas shown 
regions of string tension (Gfi) and loop size (e) that are ruled out 
by our measurement for various values of reconnection probability 
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Fig. 11. — Cosmic string constraints in the large loop case, in 
terms of string tension G/x and reconnection probability p. The 
shaded area is ruled out by our GW upper limit. 

5.3. Future Performance 

Predicting the future sensitivity of a pulsar timing 
array project to GW depends on a large number of 
poorly constrained factors. These include the level 
of pulsar-intrinsic timi ng noise in the MSP population 
jSh annon fc Co rdcs 2010), the statistical characteristics 
of the pulsar radio emission (.Oslowski et al.ii2011 ). and 
the influence of th e interstellar medium on pulsar timing 
(jColes et al.ll2010[ ). A detailed assessment of the effect 
of variou s kinds of noise on GW det ection sensitivty was 
done by iCordes fc Shannon! ()2011l ). Assuming all fac- 
tors besides experiment duration T (e.g., number of pul- 
sars, observing cadence, telescope instrumentation, etc) 
remain the same, we can roughly bound our future GW 
limits by two cases: 

• Most optimistically, the data could remain dom- 
inated primarily by "white" radiometer noise. In 
this case, limits on Ai scale in proportion to T"-3/2 
For example, Ai oc T~^^/^ for the supermassive 
black hole GWB. With an additional five years of 
data, we will reach Ai ^ 1.6 x 10^^*^ for this spec- 
trum, well into the range of expected GW ampli- 
tudes. 

• On the other hand, if the data are completely dom- 
inated by timing noise with spectral index equal to 
or steeper than the GW signal, limits improve only 
as T~^/^. Doubling the data span would only im- 
prove the result to ~ 5 x 10~^^. Note that this 
does not appear to be the situation currently for 
the majority of our pulsars. 

Intermediate cases, which are more likely to apply in 
reality than either of the two extremes above, include 
those where only some of the pulsars are timing-noise 
dominated, or have timing noise that is less red than 
the GW signal. Note that these scalings only apply 
for GW limits (that is, non-detections) - in the GW 
signal-dominated regime, and during the transition 
from limit to de t ection , different rules apply (see 
ICordes fc Shannon! 1201 If ). These estimates also do 
not take into account any future improvements in the 



experiment. In reality both telescopes have recently 
undergone a major upgrade in digital pulsar instru- 
mentation with the ins tallation of the GUP PI and 
PUPPI pulsar backends (jPemorest et al.l 120121 ). These 
provide an order of magnitude more radio bandwidth 
than was used for this analysis. NANOGrav is also 
currently in the process of expanding the number of 
pulsars monitored to ~30, as new MSPs are discovered 
by se arches including Ferm i gamma-ray source follow-up 
(e.g., 'Ransom et all 1201 If ), the Arecibo PALFA survey 
(jKaspi ct al. 20 12), and the Gree n Bank North Celestial 
Cap survey (jStovall et al.l l2012f ) . Additional improve- 
ments may come from connecting our current data set 
to historical pulsar timing data stretching back as far as 
20 years for some sources. 

6. CONCLUSIONS 

In this paper, we presented and analyzed five years 
of radio timing data on 17 pulsars taken with the two 
largest single-dish telescopes in the world. Our timing 
analysis included novel methods for dealing with time- 
variable dispersion measure and intrinsic profile shape 
evolution with frequency. We achieved sub-microsecond 
RMS timing results on all but two of our sources, and 
RMS residuals of only ~30-50 ns in the two best cases. 
We presented a new time-domain method for detecting 
and/or limiting timing flucutations of a known spectral 
shape but unknown amplitude. The key feature of this 
analysis was proper accounting for the signal power re- 
moved by the timing model fit, without requiring depen- 
dence on simulation for this step. Applied to our timing 
of J1713-I-0747, these methods set a single-pulsar limit on 
the stochastic gravitational wave background amplitude 
of ^1 < 1.1 X 10" (95%) for the expected supermassive 
black hole GWB spectrum with spectral index a = —2/3. 
We discussed how to measure cross-correlations between 
the timing of pairs of pulsars while also accounting for 
both the timing model fit and the presence of correlated 
(non-white) noise in the data. For a — —2/3, the mea- 
sured cross-correlations in our data set constrain A^ to 
(-10 ± 26) X 10"^", or alternately a 2-a upper limit of 
Al < 7.2 X 10"^^. We discussed how our measurement 
will improve with time, and suggest that prospects are 
good for obtaining astrophyiscally constraining GW lim- 
its, or possibly even a detection, over the next five years. 
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Fig. 12.— Timing summary for PSR J0030+0451, see Figure[3]for details. 
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Fig. 14. — Timing summary for PSR J1012+5307, see Figure[3]for details. 
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Timing summary for PSR J1455— 3330, see Figure |3] for details. 
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Fig. 16.— Timing summary for PSR J1600-3053, see Figure|3]for details. 
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Fig. 18.— Timing summary for PSR J1643-1224, see Figure|3]for details. 
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Fig. 19.— 



Timing summary for PSR J1744— 1134, see Figure|3]for details. 
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Fig. 20.— Timing summary for PSR J1853+1308, sec Figure[3]for details. 



CO 

■D 



CD 
CC 



CO 

■g 
en 

CD 
CC 



30 
20 
10 

-10 
-20 
-30 



3 
2 
1 

-1 
-2 
-3 




il 




2005 



2006 



2007 2008 
Year 



2009 



i T 



i if 



2005 



2006 



2007 



2008 



2009 



Year 



Year 



2010 



2010 



2010 




Fig. 21. — Timing summary for PSR B1855+09, see Figure[3]for details. 
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Fig. 22.— Timing summary for PSR J1909-3744, see Figure|3]for details. 
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Timing summary for PSR J1910-I-1256, see Figure[3]for details. 
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Fig. 25.— Timing summary for PSR B1953-I-29, see Figure[3]for details. 
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Fig. 26. — Timing summary for PSR J2145-0750, see Figure|3]for details. 
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